The present invention relates to the generation of images from projection measurements. Examples of images generated from projection measurements include two-dimensional and three-dimensional SAR (synthetic aperture radar) systems. SAR is a form of radar in which the large, highly-directional rotating antenna used by conventional radar is replaced with many low-directivity small stationary antennas scattered over some area near or around the target area. The many echo waveforms received at the different antenna positions are post-processed to resolve the target. SAR can be implemented by moving one or more antennas over relatively immobile targets, by placing multiple stationary antennas over a relatively large area, or combinations thereof. A further example of images generated from projection measurements are ISAR (inverse SAR) systems, which image objects and many features on the ground from satellites, aircraft, vehicles or any other moving platform. SAR and ISAR systems are used in detecting, locating and sometimes identifying ships, ground vehicles, mines, buried pipes, roadway faults, tunnels, leaking buried pipes, etc., as well as discovering and measuring geological features, forest features, mining volumes, etc., and general mapping. For example, as shown in FIG. 1 of U.S. Pat. No. 5,805,098 to McCorkle, hereby incorporated by reference, an aircraft mounted detector array is utilized to take ground radar measurements. Other examples of systems using projection measurements are fault inspection systems using acoustic imaging, submarine sonar for imaging underwater objects, seismic imaging system for tunnel detection, oil exploration, geological surveys, etc., and medical diagnostic tools such as sonograms, echocardiograms, x-ray CAT (computer-aided tomography) equipment and MRI (magnetic resonance imaging) equipment.
Systems which produce images from projection data generally use techniques in the time domain, where a backprojection-type algorithm is used, or frequency domain, where Fourier transforms are used. Since a Fast Fourier Transform (FFT) technique, such as a technique known as the “ω-k” implementation, requires data to be equally spaced, FFT-based techniques produce sub-optimal images when the data source is moving uncontrollably, such as an aircraft buffeted by winds or vehicles in rough terrain. Non-Uniform spacing requires a Discrete Fourier Transform (DFT) which increases computation expense relative to a backprojector technique. Also, two-dimensional FFT's are not well suited to multiprocessor-based supercomputers because they face a corner-turn interprocessor communication bottleneck.
While there are many forms of Fourier-based algorithms for SAR processing, they fall into two broad classes known as “strip-map” mode and “spot light” mode. The most robust technique is the ω-k technique, also known as seismic migration. The advantage of the ω-k algorithm over the backprojection algorithm is speed. The ω-k algorithm is an order N2 log2(N) implementation—much faster than N3 for large images and data sets.
Time domain backprojection-based techniques have been used for numerous applications, including x-ray CAT scans, MRI and sonograms. Historically, medical people have preferred backprojection because its artifact levels were lower than those using fast Fourier transform (FFT) approaches. Efforts in the past to speed up the backprojection process have focused on fast index generation. The algorithm form used by the medical industry (e.g., Star Computers) for x-ray CAT scans requires approximately 2N3 adds to form an N by N image from N projections—N3 adds for indexing operations, and N3 adds for accumulating the projections into the image. Seismologists and people using SAR have also used backprojection.
Synthetic aperture radar systems have been used in applications such as area mapping, surveillance, and target detection. The radar is usually mounted on an aircraft or a vehicle configured with transmitting and receiving antennas to transmit and measure the reflected radar signals from areas of interest. Through signal processing, the reflected radar signals along the flight path are combined to form the SAR imaging for side looking or forward looking surveillance.
SAR imaging is complex for a variety of reasons. First, the data is not inputted at equally distant (or known) points. Instead, data may be inputted in a non-uniform manner from an aircraft that is buffeted by the wind or from a ground vehicle that traverses rough ground. Therefore, motion compensation must be introduced in order to produce sharp images. Second, the subject objects need not be point sources but may be dispersive—where energy is stored and “re-radiated” over time. Ground penetrating SAR adds the complication that the media propagation velocity varies which complicates seismic processing. For many SAR applications, especially for high-resolution, ultra-wide-angle (UWA), ultra-wide-bandwidth (UWB) surveillance systems, the task is particularly problematic because the data sets are large, real-time operation is essential, and the aperture geometry is not controlled. For small aircraft buffeted by the wind can affect SAR data due to significant off-track motion and velocity changes. As a result, the data is not sampled at equally spaced intervals.
Backprojection techniques provide many advantages; including sharper images. Although prior art backprojector implementations may generate image artifacts; they are constrained to be local to the object generating the artifacts and generally lie within the theoretical sidelobes. Side lobes are the lobes of the radiation pattern that are not the main beam or lobe. In an antenna radiation pattern or beam pattern, the power density in the side lobes is generally much less than that in the main beam. It is generally desirable to minimize the sidelobe level (SLL), commonly measured in decibels relative to the peak of the main beam. The concepts of main and side lobes apply to (but are not limited to) for example, radar and optics (two specific applications of electromagnetics) and sonar. The present invention is directed to techniques which minimize the effects of theoretical sidelobes in order to provide enhanced images.
Backprojector techniques also allow for non-uniform spacing of the projection data. The non-uniform spacing is directly accounted for in the index generation, which is important when compensating for aircraft motion.
Conventional time domain image formation, or backprojection, from SAR data, is accomplished by coherently summing the sampled radar returns for each pixel. In this context, coherent summation can be thought of as time-shifting the signal obtained at each aperture position (to align them to a particular pixel) and adding across all aperture positions to integrate the value at that pixel. This time-align-and-sum sequence is repeated for every pixel in the image.
A method and system for forming images by backprojection is explained in U.S. Pat. No. 5,805,098 to McCorkle, hereby incorporated by reference as though fully rewritten herein. Specifically, FIG. 2 of the 1998 patent illustrates antennas at positions 208 along axis 204 in an array that observe pixels 202 in the near field of the array. A relative position of each pixel (q,r) with respect to each antenna position j defines a vector 206. For each pixel (q,r), the disclosed process time-shifts the signals obtained at each aperture position j (to align, or stack, them at a particular pixel location) to correct the signals for propagation time along each vector 206 and then adds across all aperture positions to integrate to the value at the pixel. Thus, signals propagating from that location are in-phase and reinforced, while signals propagating from other locations are not in phase and integrate toward zero. The image is generated by forming such a sum for each pixel as shown in equation (1A) below.
In equation (1A) below, the pixels of the image area are indexed by (q,r) and the aperture positions are indexed by j, where j=0 . . . L−1 and L is the number of elements in the aperture. If sj(t) represents the range-compensated (R2 propagation loss corrected, where R is range) voltage received at the jth aperture position as a function of time (t), zj is an aperture weighting to shape the sidelobes, for example, with a Hamming window, or to account for the aperture spacing, and Tq,r,j is the time shift necessary to align the signal received at sensor position j to the pixel at position (q,r) (a function of the round-trip time from sensor phase center to pixel position), then the value of the focused pixel at image position (q,r) is given by:
                                          f                          q              ,              r                                ⁡                      (            t            )                          =                              ∑                          j              =              0                                      L              -              1                                ⁢                                          ⁢                                    z              i                        ·                                                            s                  j                                ⁡                                  (                                      t                    +                                          T                                              q                        ,                        r                        ,                        j                                                                              )                                            .                                                          (                  1          ⁢          a                )            
Here, t describes how the focused signal at location (q,r) varies with time, and is useful for studying late-time target ringing. This description of backprojection considers the case where t is fixed for the entire image.
Accurately obtaining the time-shifted values sj(t+Tq,r,j) requires a time domain interpolation of the sampled received signals. Prior art techniques included the following steps:                1A. Up-sample and low-pass filter the received signal to produce a finer resolution signal sj.        2A. Compute the floating point indexes into the sequence s.sub.j corresponding to time t+Tq,r,j.        3A. Linearly interpolate between samples to obtain an approximation of s.sub.j (t+T.sub.q,r,j).        
The following references give an overview of the state of the art and are hereby incorporated by reference in their entireties:    1. J. McCorkle, “Focusing Of Synthetic Aperture Ultra Wideband Data,” IEEE Int'l Conf on Systems Engineering, August, 1992, p. 1-5;    2. J. McCorkle and Lam Nguyen, “Focusing of Dispersive Targets Using Synthetic Aperture Radar,” ARL-TR-305, August, 1994;    3. R. Stolt, “Migration by Fourier Transform,” Geophysics, Vol. 43, p. 23ff;    4. F. Rocca, C. Cafforio, and C. Prati, “Synthetic Aperture Radar: A New Application for Wave Equation Techniques,” Geophysical Prospecting, Vol. 37, 1989, pp. 809-30.    5. C. Cafforio, C. Prati, and F. Rocca, “SAR Data Focusing Using Seismic. Migration Techniques,” IEEE Transactions on Aerospace and Electronic Systems, Vol. AES-27, No. 2, March, 1991, pp. 194-206;    6. R. Bamler, “A Comparison of Range Doppler and Wavenumber Domain SAR. Focusing. Algorithms,” IEEE Transactions on Geoscience and Remote. Sensing, Vol. 30, No. 4, Jul. 1, 1992, pp. 706-713;    7. M. Ressler et al., “The Army Research Laboratory Ultra-Wideband Testbed Radars,” IEEE 1995 International Radar Conference, Alexandria, Va., May, 1995; and    8. L. Happ et al., “Low-Frequency Ultra-Wideband Synthetic Aperture Radar 1995 BoomSAR Tests,” IEEE 1996 National Radar Conference, Ann Arbor, Mich., May, 1996.
An example of a forward-looking Synchronous Impulse Reconstruction (SIRE) radar that can be vehicle-mounted has been designed and built by the Army Research Lab. A more complete description of the SIRE radar can be found in M. Ressler, L. Nguyen, F. Koenig, D. Wong, and G. Smith, “The Army Research Laboratory (ARL) Synchronous Impulse Reconstruction (SIRE) Forward-Looking Radar”, Proceedings of SPIE, Unmanned Systems Technology IX, April 2007, hereby incorporated by reference. The SIRE radar has two transmitters and an array of receiving antennas. The two transmitters alternatively transmit wide bandwidth impulses to illuminate the area in front of the vehicle. An array of receiving antennas measures the returned radar signals. The wide bandwidth of transmitted impulses provides the down-range resolution while the array of receiving antennas provides the cross-range resolution. It has been shown that the configuration with two transmitters located at the end of the array is the optimum configuration to achieve cross-range resolution while minimizing the number of required transmitters.
After data is acquired by the radar hardware, it is transferred to a computer for signal processing and image formation. The signal processing stages include a) self-interference extraction, b) removing radar signature distortion due to moving platform, and c) sub-band filtering. The self-interference processing step to extract the interference components from the returned radar signals and the technique to remove the phase and shape distortion in radar signals due to the motion of the radar platform are described in the publication by Lam Nguyen, entitled “Signal Processing Technique to Remove Signature Distortion in ARL Synchronous Impulse Reconstruction (SIRE) Ultra-Wideband (UWB) Radar,” Army Research Laboratory Technical Report, ARL-TR-4404, March 2008, hereby incorporated by reference.
After all the signal processing steps are applied to the returned radar signals, the processed radar range profiles may be used for forming a SAR image. In a preferred embodiment, the back-projection algorithm is utilized for the image formation step. See, John McCorkle and Lam Nguyen, “Focusing of Dispersive Targets Using Synthetic Aperture Radar,” Army Research Laboratory Report, ARL-TR-305, August 1994.
FIG. 1A illustrates an example utilizing the basic concept of the backprojection imaging algorithm. The radar is mounted on a moving platform. It transmits radar signals to illuminate the area of interest and receives return signals from the area. Using the motion of the platform, the radar collects K data records along its path (or aperture). In general the aperture could be a line, a curve, a circle, or any arbitrary shape. The receiving element k from the aperture is located at the coordinate (xR(k), yR(k), zR(k)). Far bistatic radar (the transmitting antenna is separate from the receiving antenna) the transmitting element k from the aperture is located at the coordinate (xT(k), yT(k), zT(k)). For monostatic radar (the transmitting antenna is the same as or co-located with the receiving antenna) the transmitting coordinates (xT(k), yT(k), zT(k)) would be the same as the receiving coordinates (xR(k), yR(k), zR(k)). Since the monostatic radar case is a special case of the bistatic radar configuration, the algorithm described here is applicable for both configurations. The returned radar signal at this receiving element k is sk(t). In order to form an image from the area of interest, we form an imaging grid that consists of N image pixels. Each pixel. Pi from the imaging grid is located at coordinate (xP(i), yP(i), zP(i)). The imaging grid is usually defined as a 2-D rectangular shape. In general, however, the image grid could be arbitrary. For example, a 3-D imaging grid would be formed for ground penetration radar to detect targets and structures buried underground. Another example is 3-D image of inside human body. Each measured range profile sk (t) is corrected for the R2 propagation loss, i.e. s′k(t)=R2(t)sk(t), where
      R    ⁡          (      t      )        =      ct    2  and c=2.997e8 m/sec. The backprojection value at pixel P(i) is
                                          P            ⁡                          (              i              )                                =                                    ∑                              k                =                1                            K                        ⁢                                                  ⁢                                          w                k                            ⁢                                                s                  k                  ′                                ⁡                                  (                                      f                    ⁡                                          (                                              i                        ,                        k                                            )                                                        )                                                                    ,                                  ⁢                  1          ≤          i          ≤          N                                    (        1        )            where wk is the weight factor and f(i,k) is the delay index to s′k(t) necessary to coherently integrate the value for pixel P(i) from the measured data at receiving element k.
The index is computed using the round-trip distance between the transmitting element, the image (pixel), and the receiving element. The transmitting element is located at the coordinate (xT(k), yT(k), zT(k)). The distance between the transmitting element and the image pixel P(i) is:
                                          d            1                    ⁡                      (                          i              ,              k                        )                          =                                                            [                                  (                                                                                    x                        T                                            ⁡                                              (                        k                        )                                                              -                                                                  x                        P                                            ⁡                                              (                        i                        )                                                                              )                                ]                            2                        +                                          [                                  (                                                                                    y                        T                                            ⁡                                              (                        k                        )                                                              -                                                                  y                        P                                            ⁡                                              (                        i                        )                                                                              )                                ]                            2                        +                                          [                                  (                                                                                    z                        T                                            ⁡                                              (                        k                        )                                                              -                                                                  z                        P                                            ⁡                                              (                        i                        )                                                                              )                                ]                            2                                                          (        2        )            
The distance between the receiving element and the image pixel P(i) is
                                          d            2                    ⁡                      (                          i              ,              k                        )                          =                                                            [                                  (                                                                                    x                        R                                            ⁡                                              (                        k                        )                                                              -                                                                  x                        P                                            ⁡                                              (                        i                        )                                                                              )                                ]                            2                        +                                          [                                  (                                                                                    y                        R                                            ⁡                                              (                        k                        )                                                              -                                                                  y                        P                                            ⁡                                              (                        i                        )                                                                              )                                ]                            2                        +                                          [                                  (                                                                                    z                        R                                            ⁡                                              (                        k                        )                                                              -                                                                  z                        P                                            ⁡                                              (                        i                        )                                                                              )                                ]                            2                                                          (        3        )            The total distance isd(i,k)=d1(i,k)+d2(i,k)  (4)The delay index is
                              f          ⁡                      (                          i              ,              k                        )                          =                              d            ⁡                          (                              i                ,                k                            )                                c                                    (        5        )            
FIG. 1B illustrates a typical imaging geometry for an ultra wide band forward looking (e.g., SIRE) radar. In this case, the radar is configured in forward-looking mode instead of side-looking mode as illustrated in FIG. 1A. In this forward-looking mode, the radar travels and radiates energy in the same direction. The general backprojection algorithm described applies to the embodiment of FIG. 1B. As seen in FIG. 1B, the radar travels in parallel to the x-axis. The backprojection image formation is combined with the mosaic technique. The large area image is divided into sub-images. The size of each sub-image may be, for example, 25 m in cross-range and only 2 m in down-range (x-axis direction). The radar starts at coordinate A, which is 20 m from sub-image 1, and illuminates the entire image area to the right.
The following is a description of the platform 10 in FIG. 1B as it passes four sequential positions 10A, 10B 10C & 10D located at x-coordinates A, B, C & D, respectively. The formation of the first sub-image begins when platform 10 is at the coordinate A, 20 meters from the block labeled “1st sub-image.” As platform 10 travels in the x direction (as shown in FIG. 1B), signals emitted from platform 10 illuminates an entire image area to the right of platform 10, and the reflected signals are received by an array of 16 physical receiving antennas 11 positioned on the front of the platform 10. Formation of the first sub-image ends when platform 10 reaches coordinate C, at approximately 8 m from the block labeled “1st sub-image.” Accordingly, the radar signal data for the first (full-resolution) sub-image is received as radar platform 10 travels a distance of 12 meters (20 m−8 m=12 m) from coordinates A to C, for formation of a two dimensional (2D) aperture.
The distance traveled during the formation of the two-dimensional (2-D) aperture is represented by an arrow in FIG. 1B labeled “Aperture 1.” When the platform 10 reaches coordinate B, a distance of 2 meters from coordinate A in FIG. 1B, the formation of the “2nd sub-image” begins, and as the platform 10 travels to coordinate D, it uses the received data to form a second 2-D aperture. The distance traveled by platform 10 is represented by an arrow in FIG. 1B labeled “Aperture 2.”Note that the two apertures are overlapped by 10 m and the second aperture is “advanced” by 2 m with respect to the first aperture. Sub-images 1 and 2 are formed from the 2-D apertures using the same length of travel (12 meters) of the radar. This process is applied to ensure that image pixels have almost the same (within a specified tolerance) resolution across the entire large area. The sub-images are formed from the radar range profiles using the back-projection algorithm.
FIG. 2 shows the back-projection algorithm applied to form a sub-image. The procedure mathematically described with respect to FIG. 1A in the above paragraphs may also be applied to this imaging scenario. In this case, the radar aperture is a rectangular array that is formed by an array of 16 receiving elements (that spans 2 meters) and the forward motion of the platform (12 meter for forming each sub-image). The imaging grid in this case is defined as a rectangular array of 25×2 meter.
FIG. 3 shows a SAR image formed using the above algorithm using simulated data of two targets (points). The image, is displayed using 40 dB of dynamic range. However, “energy” from the two point targets is spread throughout the image and creates severe sidelobes. There are two sources that generate the imaging artifacts in this case: a) aperture aliasing (small aperture compared to the large image cross-range swath), and b) the errors from the position measurements system. In reality, there are many other sources that contribute to the noise floor of the resulting image. This created a challenging problem for the detection of targets of smaller amplitudes since they might be obscured or even embedded in the noise floor.
The term “noise” as used herein relates to image noise. There are many sources that cause noise in the resulting image. Noise can be divided into two categories: additive noise and multiplicative noise. System noise, thermal noise, quantization noise, self-interference noise, radio frequency interference (RFI) noise are some examples of the additive noise. Multiplicative noise is much more difficult to deal with since it is data dependent. Some sources that cause multiplicative noise include: timing jitter in data sampling, small aperture size compared to image area, the under-sampling of aperture samples, the non-uniform spacing between aperture samples, errors in position measurement system, etc. Multiplicative noise results in undesired sidelobes that create high noise floor in the image and thus limit the ability to detect targets with smaller amplitudes.
Radar and other imaging systems currently suffer various noise sources that prevent the generation of very high contrast images. As a result, difficult targets (with low amplitudes) are often obscured or even embedded, in the noise level of the image background. Moreover, sidelobes from large targets are mistaken as targets of interest. Recently the ARL has designed and built a new ultra-wideband imaging system for the detection of difficult targets. Currently, there exists a need for an improved signal processing technique which reduces unwanted noise and enhances image reproduction.
Stepped-Frequency radar technology comprises hardware and software subsystems that have attained advanced levels of technological maturity. The popularity of the paradigm attests to both its practicality and its analytic tractability. Stepped-frequency radar systems achieve high resolution in the time domain by effectively measuring the frequency-domain response over a large bandwidth and then transforming from frequency to time. In particular, a series of pulses is transmitted, with each pulse centered at a specific frequency. The reflected waves from a target are coherently sampled, and since the transmitted frequencies are evenly-spaced, a high resolution time-domain, waveform can readily be constructed via a Fast Fourier Transform (FFT). The time-domain resolution of the system is determined by the total bandwidth (i.e. the number of frequency steps and the frequency step size). The duration of the high-resolution signature will be fixed by the pulse width and the frequency step size.
When operating at lower frequencies, it is possible for stepped-frequency radars to interfere with indigenous radio frequency (RF) transmissions, such as television stations, cell phones, and air traffic control. Hence, it becomes necessary for a radar operator to “notch-out” certain transmit frequencies, thereby producing unwanted artifacts in the final radar signature. These effects are well-documented in the literature, and efforts have been undertaken to mitigate them, with varying degrees of success. 2. STEPPED FREQUENCY RECURSIVE SIDELOBE MINIMIZATION
Non-linear techniques for reducing sidelobe artifacts have existed for several years. For example, non-linear apodization techniques, developed by Stankwitz et al. 6-7, exploit the time- and frequency-domain characteristics of specific weighting windows to reduce sidelobe artifacts in synthetic aperture radar (SAR) imagery. Another recently presented technique, recursive sidelobe minimization (RSM), operates strictly in the time-domain, and it combines concepts from the theories of both apodization and compressive sensing 8-10 to reduce image artifacts. Unfortunately, this restriction of the current RSM formulation to the time domain precludes elimination of artifacts introduced by frequency notching.